Demonstration for Comparing the MSIS model versions in GEODYN

  • This notebook shows the output for running GEODYN with the three latest versions of the MSIS series of upper atmospheric models: MSIS_86, MSIS_00, and MSIS_2.0.

  • The run parameters for the GEODYN run demoed here are: - Satellite:             Starlette          - Observation Datatype:  SLR          - Arc time:              14 days (2003/09/14 - 2003/09/28)          - Gen. Accelerations:    Off          - Density Model:         NRLMSISe86, NRLMSISe00, MSIS 2.0

  • The challenges that were overcome to include these three models into versions of GEODYN are detailed in an attached document.

Input Starlette-SLR run parameters:

  • input the parameters that make up the starlette-SLR run type in GEODYN

[1]:
%load_ext autoreload
%autoreload 2

###################################
                                  #
# CHOOSE Models to compare:       #
m1 = 'msis86'                     #
m2 = 'msis00'                     #
m3 = 'msis2'                      #
###################################



##################################
# INPUT PARAMETERS:
##################################
sat = 'st'
arc = '030914_2wk'
grav_id ='goco05s'
local_path = '/data/geodyn_proj/analysis/starlette_analysis/'
Accel_Status = 'acceloff'
SAT_ID = 7501001



##################################
# PATH TO DENSITY MODEL RUNS:
##################################

msis86_model = 'msis86'
path_to_msis86 = '/data/geodyn_proj/runs_geodyn'+'/results/' + sat+'/' +msis86_model+'/'+  msis86_model+'_'+ Accel_Status +'/'

msis2_model = 'msis2'
path_to_msis2 = '/data/geodyn_proj/runs_geodyn'+'/results/' + sat+'/' +msis2_model+'/'+  msis2_model+'_'+ Accel_Status +'/'

msis00_model = 'msis00'
path_to_msis00 = '/data/geodyn_proj/runs_geodyn'+'/results/' + sat+'/' +msis00_model+'/'+  msis00_model+'_'+ Accel_Status +'/'


if m1 == 'msis86':
    path_to_m1 = path_to_msis86
elif m1 == 'msis00':
    path_to_m1 = path_to_msis00
elif m1 == 'msis2':
    path_to_m1 = path_to_msis2


if m2 == 'msis86':
    path_to_m2 = path_to_msis86
elif m2 == 'msis00':
    path_to_m2 = path_to_msis00
elif m2 == 'msis2':
    path_to_m2 = path_to_msis2

if m3 == 'msis86':
    path_to_m3 = path_to_msis86
elif m3 == 'msis00':
    path_to_m3 = path_to_msis00
elif m3 == 'msis2':
    path_to_m3 = path_to_msis2



import sys
sys.path.insert(0, '/data/geodyn_proj/analysis/util_funcs/py_read_geodyn_output/')
from a_ReadGEODYN import Read_GEODYN_func



(AdjustedParams_m1,
 Trajectory_m1,
 Density_m1,
 ResidsObs_m1,
 ResidsSummStation_m1,
 ResidsMeasSumm_m1,
 Stats_m1)  = Read_GEODYN_func(arc,
                            sat,
                            SAT_ID,
                            grav_id,
                            m1,
                            local_path,
                            path_to_m1,
                            False,
                            2003,
                            'SLR',
                            verbose_loading=False,
                            Verbose_Stats=True)

(AdjustedParams_m2,
 Trajectory_m2,
 Density_m2,
 ResidsObs_m2,
 ResidsSummStation_m2,
 ResidsMeasSumm_m2,
 Stats_m2)  = Read_GEODYN_func(arc,
                            sat,
                            SAT_ID,
                            grav_id,
                            m2,
                            local_path,
                            path_to_m2,
                            False,
                            2003,
                            'SLR',
                            verbose_loading=False,
                            Verbose_Stats=True)


(AdjustedParams_m3,
 Trajectory_m3,
 Density_m3,
 ResidsObs_m3,
 ResidsSummStation_m3,
 ResidsMeasSumm_m3,
 Stats_m3)  = Read_GEODYN_func(arc,
                            sat,
                            SAT_ID,
                            grav_id,
                            m3,
                            local_path,
                            path_to_m3,
                            False,
                            2003,
                            'SLR',
                            verbose_loading=False,
                            Verbose_Stats=True)

              Loading  /data/geodyn_proj/runs_geodyn/results/st/msis86/msis86_acceloff/


+============================== Msis86 Run Details ==============================+
     Density model: msis86
     Satellite: starlette
     Data type: SLR

 CONVERGENCE WITHIN  2.0 PERCENT AFTER  6 ITERATIONS

 COORDINATE SYSTEM REFERENCE EPOCH 1030914    0  0.0000

 START 1030914    0  0.0000EPOCH 1030914    0  0.0000  END 1030928    0  0.0000
 INTEGRATION STEP SIZE -- ORBIT(SEC) =     15.000  VAR.EQ.(SEC) =     15.000

 SAT. ID = 7501001  AREA(M**2) =     0.04524  MASS(KG) =      47.250
    81286 ORBIT INTEGRATION STEPS  --  AVERAGE NUMBER OF CORRECTIONS = 1.0000
    81286 VARIATIONAL EQUATION INTEGRATION STEPS --  ORBIT STARTS IN SUNLIGHT
      390 TRANSITIONS BETWEEN SHADOW AND SUNLIGHT  --  28.0 PERCENT IN SHADOW
 X POS  =  -1228611.047711857     M      S.M.A. =   7329875.737190500     M
 Y POS  =  -4721160.098619740     M      ECCEN. =  0.2046448692470624E-01
 Z POS  =   5602824.670083221     M      INCLIN =   49.80801875659012     DEG
 X VEL  =   6644.663061993014     M/S    NODE   =   151.4090092660359     DEG
 Y VEL  =  -2852.907078978589     M/S    PERG   =   326.8920903447495     DEG
 Z VEL  =  -798.7567429058066     M/S    MEAN   =   130.4976391805683     DEG
 RMS POS=            0.051571     M      RMS VEL=            0.000050     M/S

 THE FOLLOWING ARE GEOCENTRIC LATITUDE AND GEOCENTRIC LONGITUDE AND AN
 APPROXIMATE HEIGHT FROM THE SEMI-MAJOR AXIS OF THE CENTRAL BODY
  LATITUDE =  48.954 DEG   LONGITUDE = 262.855 DEG   HEIGHT =         1050893.385 M

 APOGEE =         1110727.926 METERS     ASC NODE RATE =    -3.957719 DEG/DAY
 PERIGEE=          810723.634 METERS     ARG PERG RATE =     3.315632 DEG/DAY
 PERIOD =             104.089 MINUTES    PERIOD   RATE =    -0.000002 SEC/DAY
 DRAG   =               0.055 M/DAY**2   S-M AXIS RATE =    -0.001199 M/DAY
+==============================================================================+
              Loading  /data/geodyn_proj/runs_geodyn/results/st/msis00/msis00_acceloff/


+============================== Msis00 Run Details ==============================+
     Density model: msis00
     Satellite: starlette
     Data type: SLR

 CONVERGENCE WITHIN  2.0 PERCENT AFTER  6 ITERATIONS

 COORDINATE SYSTEM REFERENCE EPOCH 1030914    0  0.0000

 START 1030914    0  0.0000EPOCH 1030914    0  0.0000  END 1030928    0  0.0000
 INTEGRATION STEP SIZE -- ORBIT(SEC) =     15.000  VAR.EQ.(SEC) =     15.000

 SAT. ID = 7501001  AREA(M**2) =     0.04524  MASS(KG) =      47.250
    81286 ORBIT INTEGRATION STEPS  --  AVERAGE NUMBER OF CORRECTIONS = 1.0000
    81286 VARIATIONAL EQUATION INTEGRATION STEPS --  ORBIT STARTS IN SUNLIGHT
      390 TRANSITIONS BETWEEN SHADOW AND SUNLIGHT  --  28.0 PERCENT IN SHADOW
 X POS  =  -1228611.000644256     M      S.M.A. =   7329875.737425419     M
 Y POS  =  -4721160.116960314     M      ECCEN. =  0.2046448684995622E-01
 Z POS  =   5602824.663907633     M      INCLIN =   49.80801876199792     DEG
 X VEL  =   6644.663072193405     M/S    NODE   =   151.4090092804378     DEG
 Y VEL  =  -2852.907047211658     M/S    PERG   =   326.8920910634153     DEG
 Z VEL  =  -798.7567797112002     M/S    MEAN   =   130.4976388422482     DEG
 RMS POS=            0.050316     M      RMS VEL=            0.000049     M/S

 THE FOLLOWING ARE GEOCENTRIC LATITUDE AND GEOCENTRIC LONGITUDE AND AN
 APPROXIMATE HEIGHT FROM THE SEMI-MAJOR AXIS OF THE CENTRAL BODY
  LATITUDE =  48.954 DEG   LONGITUDE = 262.855 DEG   HEIGHT =         1050893.385 M

 APOGEE =         1110727.926 METERS     ASC NODE RATE =    -3.957719 DEG/DAY
 PERIGEE=          810723.635 METERS     ARG PERG RATE =     3.315632 DEG/DAY
 PERIOD =             104.089 MINUTES    PERIOD   RATE =    -0.000001 SEC/DAY
 DRAG   =               0.050 M/DAY**2   S-M AXIS RATE =    -0.001103 M/DAY
+==============================================================================+
              Loading  /data/geodyn_proj/runs_geodyn/results/st/msis2/msis2_acceloff/


+============================== Msis2 Run Details ==============================+
     Density model: msis2
     Satellite: starlette
     Data type: SLR

 CONVERGENCE WITHIN  2.0 PERCENT AFTER  6 ITERATIONS

 COORDINATE SYSTEM REFERENCE EPOCH 1030914    0  0.0000

 START 1030914    0  0.0000EPOCH 1030914    0  0.0000  END 1030928    0  0.0000
 INTEGRATION STEP SIZE -- ORBIT(SEC) =     15.000  VAR.EQ.(SEC) =     15.000

 SAT. ID = 7501001  AREA(M**2) =     0.04524  MASS(KG) =      47.250
    81286 ORBIT INTEGRATION STEPS  --  AVERAGE NUMBER OF CORRECTIONS = 1.0000
    81286 VARIATIONAL EQUATION INTEGRATION STEPS --  ORBIT STARTS IN SUNLIGHT
      390 TRANSITIONS BETWEEN SHADOW AND SUNLIGHT  --  28.0 PERCENT IN SHADOW
 X POS  =  -1228611.002860558     M      S.M.A. =   7329875.737429938     M
 Y POS  =  -4721160.115078840     M      ECCEN. =  0.2046448656421785E-01
 Z POS  =   5602824.663322780     M      INCLIN =   49.80801876249247     DEG
 X VEL  =   6644.663073195048     M/S    NODE   =   151.4090092811391     DEG
 Y VEL  =  -2852.907048276599     M/S    PERG   =   326.8920909993722     DEG
 Z VEL  =  -798.7567790814259     M/S    MEAN   =   130.4976389110996     DEG
 RMS POS=            0.050103     M      RMS VEL=            0.000049     M/S

 THE FOLLOWING ARE GEOCENTRIC LATITUDE AND GEOCENTRIC LONGITUDE AND AN
 APPROXIMATE HEIGHT FROM THE SEMI-MAJOR AXIS OF THE CENTRAL BODY
  LATITUDE =  48.954 DEG   LONGITUDE = 262.855 DEG   HEIGHT =         1050893.383 M

 APOGEE =         1110727.924 METERS     ASC NODE RATE =    -3.957719 DEG/DAY
 PERIGEE=          810723.637 METERS     ARG PERG RATE =     3.315632 DEG/DAY
 PERIOD =             104.089 MINUTES    PERIOD   RATE =    -0.000001 SEC/DAY
 DRAG   =               0.050 M/DAY**2   S-M AXIS RATE =    -0.001106 M/DAY
+==============================================================================+

Coefficient of Drag Adjustments

Details and Background:

  • the \(C_d\) is an adjusted parameter in GEODYN. Each iteration, it is estimated and adjusted to improve the OD simulation

  • \(C_d\) is treated as a constant in GEODYN (for any iteration). The factor \(C_d\) varies slightly with satellite shape and atmospheric composition. However, for any geodetically useful satellite, it may be treated as a satellite dependent constant. If the Cd is time dependent, then it is constant over the designated time intervals.

  • $ \bar`{A}_D = -:nbsphinx-math:frac{1}{2}` C_d \frac{A_s}{m_s} \rho_d v_r :nbsphinx-math:`bar`{v}_r $

    • \(A_D\) acceleration due to atmospheric drag force

    • \(C_d\) is the satellite drag coefficient

    • \(A_s\) is the cross-sectional area of the satellite

    • \(m_s\) is the mass of the satellite

    • \(\rho_d\)is the-density of the atmosphere at the satellite position

    • \(\bar{v}_r\) is the velocity vector of the satellite relative to the atmosphere, and

    • \(v_r\) is the magnitude of the velocity vector, vr .

Adjusted Cd Plot (timeseries)

[2]:
labels = list(AdjustedParams_m1[5][SAT_ID]['0CD'].keys())
val_list = []
for i in AdjustedParams_m1[5][SAT_ID]['0CD'].keys():
    val_list.append(AdjustedParams_m1[5][SAT_ID]['0CD'][i]['CURRENT_VALUE'])



import plotly.graph_objects as go
import numpy as np
from plotly.offline import plot, iplot
#     %matplotlib inline
from plotly.subplots import make_subplots


last_iter = 5
which_stat = 'CURRENT_VALUE' # 2 is the current val

labels = list(AdjustedParams_m1[5][SAT_ID]['0CD'].keys())


fig = make_subplots(
    rows=1, cols=1,
#         subplot_titles=("Time Dependent Drag Coefficient ",
                    )

val_list = []
for i in AdjustedParams_m1[last_iter][SAT_ID]['0CD'].keys():
    val_list.append(AdjustedParams_m1[last_iter][SAT_ID]['0CD'][i][which_stat])
fig.add_trace(go.Scatter(x=labels,
                                 y=val_list,
                                  name= m1,
                               mode='markers',
                                    marker=dict(
                                    size=10,),
                                ), row=1, col=1,
                                   )


val_list = []
for i in AdjustedParams_m2[last_iter][SAT_ID]['0CD'].keys():
    val_list.append(AdjustedParams_m2[last_iter][SAT_ID]['0CD'][i][which_stat])
fig.add_trace(go.Scatter(x=labels,
                                 y=val_list,
                                  name= m2,
                               mode='markers',
                                    marker=dict(
                                    size=8,),
                                ), row=1, col=1,
                                   )

val_list = []
for i in AdjustedParams_m3[last_iter][SAT_ID]['0CD'].keys():
    val_list.append(AdjustedParams_m3[last_iter][SAT_ID]['0CD'][i][which_stat])
fig.add_trace(go.Scatter(x=labels,
                                 y=val_list,
                                  name= m3,
                               mode='markers',
                                    marker=dict(
                                    size=8,),
                                ), row=1, col=1,
                                   )

fig.update_yaxes( title="Cd (Unitless)",exponentformat= 'power',row=1, col=1)
fig.update_xaxes( title="Date", row=1, col=1)

fig.update_layout(title="Time Dependent Drag Coefficient ")
fig.update_layout(
    font=dict(
#             family="Courier New, monospace",
        size=14,
             ),)
fig.update_layout(
    autosize=False,
    width=600,
    height=500,
)
iplot(fig)

Cd Analysis

A few key features can be seen in the above plot: 1. When there is some sort of storm event early in the arc, the Cd responds by correcting the mean density provided by the different models. 2. The MSISe00 corrections are less at this set of altitudes as opposed to those of m86 and m2.0. 3. Overall, though from a POD perspective, they seem to do a good job accounting for density variations. - Frank states that the onus is to show there is some real storm event going on over these several days.

TODO: add storm event from the Kp/Ap and F10.7 timeseries.

Density Output

  • The density output is sourced from the DRAG.f90 subroutine.

[3]:
import plotly.graph_objects as go
import numpy as np
from plotly.offline import plot, iplot
#     %matplotlib inline
from plotly.subplots import make_subplots
import plotly.express as px

col1 = px.colors.qualitative.Plotly[0]
col2 = px.colors.qualitative.Plotly[1]
col3 = px.colors.qualitative.Plotly[2]

data_nums_1 = 1000
data_nums_2 = 100

fig = make_subplots(rows=2, cols=1,
    subplot_titles=("First 4 Minutes","Every 100th point"),
                   )

fig.add_trace(go.Scatter(x=Density_m1['Date'][:data_nums_1],
                                 y=Density_m1['rho (kg/m**3)'][:data_nums_1],
                                    name=m1,
                                     mode='markers',
                                    marker=dict(
                                    size=4,
                                    ),
                                   ),
                                    row=1, col=1,
                                   )


fig.add_trace(go.Scatter(x=Density_m2['Date'][:data_nums_1],
                         y=Density_m2['rho (kg/m**3)'][:data_nums_1],
                         name= m2,
                         mode='markers',
                         marker=dict(
                         size=2,),
                         ),
                        row=1, col=1,
                         )

fig.add_trace(go.Scatter(x=Density_m3['Date'][:data_nums_1],
                         y=Density_m3['rho (kg/m**3)'][:data_nums_1],
                         name= m3,
                         mode='markers',
                         marker=dict(
                         size=2,),
                         ),
                          row=1, col=1,
                         )
fig.add_trace(go.Scatter(x=Density_m1['Date'][::data_nums_2],
                                 y=Density_m1['rho (kg/m**3)'][::data_nums_2],
                                 name = m1,
                                 mode='markers',
                                    marker=dict(color=col1,
                                    size=3,
                                    ),
                                  showlegend=False,

                                   ),
                                   row=2, col=1,
                                   )

fig.add_trace(go.Scatter(x=Density_m2['Date'][::data_nums_2],
                         y=Density_m2['rho (kg/m**3)'][::data_nums_2],
                         name= m2,
                         mode='markers',
                         marker=dict(color=col2,
                         size=2,),
                                  showlegend=False,

                         ),
                        row=2, col=1,
                         )

fig.add_trace(go.Scatter(x=Density_m3['Date'][::data_nums_2],
                         y=Density_m3['rho (kg/m**3)'][::data_nums_2],
                         name= m3,
                         mode='markers',
                         marker=dict(color=col3,
                         size=2,),
                                  showlegend=False,

                         ),
                          row=2, col=1,
                         )

fig.update_layout(
    title="Density along Starlette Orbit",
    )

fig.update_yaxes(type="log", exponentformat= 'power',)
fig.update_layout(legend= {'itemsizing': 'constant'})

fig.update_yaxes(title_text="kg/m^3", row=1, col=1)
fig.update_yaxes(title_text="kg/m^3", row=2, col=1)

fig.update_xaxes(title_text="Date", row=1, col=1)
fig.update_xaxes(title_text="Date", row=2, col=1)

fig.update_layout(legend= {'itemsizing': 'constant'})
#     fig.update_xaxes(tickangle=45)


fig.update_layout(
    autosize=False,
    width=900,
    height=700,
)
fig.update_layout(
    font=dict(
#             family="Courier New, monospace",
        size=18,
             ),)

iplot(fig)




Density output Analysis

  • Starlette is in somewhat of an eccentric orbit, roughly 800 km x 1100 km orbit - so that accounts for the change in density along the orbit.

  • while this is helpful to see, one really needs to demonstrate these results over multiple arcs.

  • The relative densities are different at Apogee and perigee.

    • Apogee : MSISe86 and MSISe00 have the higher densities

    • Perigee : MSIS2.0 has the higher density

TODO: run this over multiple arcs and add an orbit averaged plot

[ ]:

Residuals – obervations

Details:

  • the residuals being output by geodyn should be rescaled to CM instead of METER. Most of your stations have cm level (or better) data precision.

  • The mean in the residuals (or the mean by station) is not important - it is very small in these runs - unless you turn out to have a station with a large systematic range bias - due to some temporary anomaly.

    • If you look at the setups, you’ll see sets of “MBIAS” cards – these adjust different sets of range biases by station per arc or by station per data pass based on recommendations from the ILRS. They are adjusted as a normal part of the POD process.

[4]:
# Rescale residuals to CM instead of M:

ResidsObs_m1['Residual'] = ResidsObs_m1['Residual']*(1e2)
ResidsObs_m2['Residual'] = ResidsObs_m2['Residual']*(1e2)
ResidsObs_m3['Residual'] = ResidsObs_m3['Residual']*(1e2)
[5]:
import plotly.graph_objects as go
import numpy as np
from plotly.offline import plot, iplot
#     %matplotlib inline
import pandas as pd

# stations_unique = []
# for station in ResidsObs_m3['track_1']:
#     if station not in stations_unique:
#         stations_unique.append(station)

# elev_angle = 40

fig = go.Figure()
# for station in stations_unique:
#     index_station = ResidsObs_m3.track_1==station
#     index_elev = ResidsObs_m3.Elev1>=elev_angle

fig.add_trace(go.Scatter(x=ResidsObs_m1['Date'],
                         y=ResidsObs_m1['Residual'].values.astype(float),
                         name= m1,
                         mode='markers',
                         marker=dict(
                         size=3,),
                         ),
                         )
fig.add_trace(go.Scatter(x=ResidsObs_m2['Date'],
                         y=ResidsObs_m2['Residual'].values.astype(float),
                         name= m2,
                         mode='markers',
                         marker=dict(
                         size=3,),
                         ),
                         )
fig.add_trace(go.Scatter(x=ResidsObs_m3['Date'],
                         y=ResidsObs_m3['Residual'].values.astype(float),
                         name= m3,
                         mode='markers',
                         marker=dict(
                         size=3,),
                         ),
                         )

fig.update_layout(
    title='Observation Residuals',
    yaxis_title="Residuals (cm)",
    xaxis_title="Date",
    )
fig.update_layout(legend= {'itemsizing': 'constant'})

# fig.update_yaxes(type="log", exponentformat= 'power')
fig.update_layout(
    font=dict(
#             family="Courier New, monospace",
        size=17,
             ),)

iplot(fig)

Analysis

  • You see the bowtie effect in the residuals above. That is a standard effect of least squares - where fits seem the best in the center of the arc.

Residuals – measurement summary

[6]:
# Rescale to CM instead of M
ResidsMeasSumm_m1['MEAN'] = ResidsMeasSumm_m1['MEAN']*1e2
ResidsMeasSumm_m2['MEAN'] = ResidsMeasSumm_m2['MEAN']*1e2
ResidsMeasSumm_m3['MEAN'] = ResidsMeasSumm_m3['MEAN']*1e2

ResidsMeasSumm_m1['WTD-MEAN'] = ResidsMeasSumm_m1['WTD-MEAN']*1e2
ResidsMeasSumm_m2['WTD-MEAN'] = ResidsMeasSumm_m2['WTD-MEAN']*1e2
ResidsMeasSumm_m3['WTD-MEAN'] = ResidsMeasSumm_m3['WTD-MEAN']*1e2
[7]:
import plotly.graph_objects as go
from plotly.subplots import make_subplots
import plotly.express as px

col1 = px.colors.qualitative.Plotly[0]
col2 = px.colors.qualitative.Plotly[1]
col3 = px.colors.qualitative.Plotly[2]

mark_size = 10
fig = make_subplots(rows=2, cols=1,
    subplot_titles=("Mean","RMS", ),
                   )
fig.add_trace(go.Scatter(x=ResidsMeasSumm_m1['Iter'][1:].astype(int),
                         y=ResidsMeasSumm_m1['MEAN'][1:].values.astype(float),
                         name= m1,
                         mode='markers',
                         marker=dict(color = col1,
                         size=mark_size,),
                         ),
                          row=1, col=1,
                         )
fig.add_trace(go.Scatter(x=ResidsMeasSumm_m2['Iter'][1:].astype(int),
                         y=ResidsMeasSumm_m2['MEAN'][1:].values.astype(float),
                         name= m2,
                         mode='markers',
                         marker=dict(color = col2,
                         size=mark_size,),
                         ),
                         row=1, col=1,
                         )
fig.add_trace(go.Scatter(x=ResidsMeasSumm_m3['Iter'][1:].astype(int),
                         y=ResidsMeasSumm_m3['MEAN'][1:].values.astype(float),
                         name= m3,
                         mode='markers',
                         marker=dict(color = col3,
                         size=mark_size,),
                         ),
                         row=1, col=1,
                         )

fig.add_trace(go.Scatter(x=ResidsMeasSumm_m1['Iter'][1:].astype(int),
                         y=ResidsMeasSumm_m1['RMS'][1:].values.astype(float),
                         name= m1,
                         mode='markers',
                         marker=dict(color = col1,
                         size=mark_size,),
                         showlegend=False,
                         ),
                          row=2, col=1,
                         )
fig.add_trace(go.Scatter(x=ResidsMeasSumm_m2['Iter'][1:].astype(int),
                         y=ResidsMeasSumm_m2['RMS'][1:].values.astype(float),
                         name= m2,
                         mode='markers',
                         marker=dict(color = col2,
                         size=mark_size,),
                         showlegend=False,
                         ),
                         row=2, col=1,
                         )
fig.add_trace(go.Scatter(x=ResidsMeasSumm_m3['Iter'][1:].astype(int),
                         y=ResidsMeasSumm_m3['RMS'][1:].values.astype(float),
                         name= m3,
                         mode='markers',
                         marker=dict(color = col3,
                         size=mark_size,),
                         showlegend=False,
                         ),
                         row=2, col=1,
                         )
fig.update_layout(
    title="Residual Measurement Summary",
#     yaxis_title="Root mean Square",
#     xaxis_title="Iteration",
    )
fig.update_yaxes(title_text="Residual (cm)", row=1, col=1)
fig.update_yaxes(title_text="RMS", row=2, col=1)

fig.update_xaxes(title_text="Iterations", row=1, col=1)
fig.update_xaxes(title_text="Iterations", row=2, col=1)

fig.update_layout(legend= {'itemsizing': 'constant'})
fig.update_xaxes(tickangle=0)

fig.update_layout(
    autosize=False,
    width=650,
    height=700,
)
fig.update_layout(
    font=dict(
#             family="Courier New, monospace",
        size=15,
             ),)

# fig.update_layout(showlegend=False)

fig.show()

Residuals – station summary

  • I have been told that these are not helpful for assessing density models.

  • unless to identify a really bad station (they can crop up from time to time because of a station anomaly). If one finds a bad station, then one can just delete the data. In this example, it looks like Riyadh (RIYA7832) seems to have a mean bias of 2 cm?. That is a bit large and might be a sign of station coordinate error or some other station-dependent error. One could solve for an MBIAS for that station over the arc to take care of that small problem - you would have to check that an MBIAS for that station is not already adjusted.

[8]:
# import plotly.graph_objects as go
# import numpy as np
# from plotly.offline import plot, iplot
# from plotly.subplots import make_subplots
# import pandas as pd
# import plotly.express as px


# col1 = px.colors.qualitative.Plotly[0]
# col2 = px.colors.qualitative.Plotly[1]
# col3 = px.colors.qualitative.Plotly[2]

# # mark_size = 10
# fig = make_subplots(rows=2, cols=1,
#     subplot_titles=("Mean","RMS", ),
#                   vertical_spacing=0.29 , )

# select_iter = ResidsSummStation_m3.Iter == 5
# # fig = go.Figure()
# mark_size = 10


# fig.add_trace(go.Scatter(x=ResidsSummStation_m1['STATION'][select_iter],
#                          y=ResidsSummStation_m1['MEAN'][select_iter].values.astype(float),
#                          name= m1,
#                          mode='markers',
#                          marker=dict(color =col1,
#                          size=mark_size,),
#                          ),
#                         row=1, col=1,
#                          )
# fig.add_trace(go.Scatter(x=ResidsSummStation_m2['STATION'][select_iter],
#                          y=ResidsSummStation_m2['MEAN'][select_iter].values.astype(float),
#                          name= m2,
#                          mode='markers',
#                          marker=dict(color =col2,
#                          size=mark_size,),
#                          ),
#                         row=1, col=1,
#                          )
# fig.add_trace(go.Scatter(x=ResidsSummStation_m3['STATION'][select_iter],
#                          y=ResidsSummStation_m3['MEAN'][select_iter].values.astype(float),
#                          name= m3,
#                          mode='markers',
#                          marker=dict(color =col3,
#                          size=mark_size,),
#                          ),
#                         row=1, col=1,
#                          )


# fig.add_trace(go.Scatter(x=ResidsSummStation_m1['STATION'][select_iter],
#                          y=ResidsSummStation_m1['RMS'][select_iter].values.astype(float),
#                          name= m1,
#                          mode='markers',
#                          marker=dict(color=col1,
#                          size=mark_size,),
#                          showlegend=False,
#                          ),
#                         row=2, col=1,
#                          )
# fig.add_trace(go.Scatter(x=ResidsSummStation_m2['STATION'][select_iter],
#                          y=ResidsSummStation_m2['RMS'][select_iter].values.astype(float),
#                          name= m2,
#                          mode='markers',
#                          marker=dict(color=col2,
#                          size=mark_size,),
#                          showlegend=False,
#                          ),
#                         row=2, col=1,
#                          )
# fig.add_trace(go.Scatter(x=ResidsSummStation_m3['STATION'][select_iter],
#                          y=ResidsSummStation_m3['RMS'][select_iter].values.astype(float),
#                          name= m3,
#                          mode='markers',
#                          marker=dict(color=col3,
#                          size=mark_size,),
#                          showlegend=False,
#                          ),
#                         row=2, col=1,
#                          )

# fig.update_layout(
#     title="Stats -- Residual Summary by Station",
# #     yaxis_title="Root mean Square",
# #     xaxis_title="Stations",
#     )

# fig.update_yaxes(title_text="mean residuals", row=1, col=1)
# fig.update_yaxes(title_text="rms", row=2, col=1)

# fig.update_xaxes(title_text="Stations", row=1, col=1)
# fig.update_xaxes(title_text="Stations", row=2, col=1)

# fig.update_layout(legend= {'itemsizing': 'constant'})
# fig.update_xaxes(tickangle=45)

# fig.update_yaxes(type="log", exponentformat= 'power')
# fig.update_layout(
#     autosize=False,
#     width=900,
#     height=700,
# )
# fig.update_layout(
#     font=dict(
# #             family="Courier New, monospace",
#         size=17,
#              ),)

# iplot(fig)


# import plotly.graph_objects as go
# import numpy as np
# from plotly.offline import plot, iplot
# # %matplotlib inline
# import pandas as pd
# from plotly.subplots import make_subplots
# import plotly.express as px


# col1 = px.colors.qualitative.Plotly[0]
# col2 = px.colors.qualitative.Plotly[1]
# col3 = px.colors.qualitative.Plotly[2]

# # mark_size = 10
# fig = make_subplots(rows=2, cols=1,
#     subplot_titles=("Weighted Mean","Weighted RMS"),
#                   vertical_spacing=0.29 ,
#                    )

# select_iter = ResidsSummStation_m3.Iter == 5
# # fig = go.Figure()
# mark_size = 10


# fig.add_trace(go.Scatter(x=ResidsSummStation_m1['STATION'][select_iter],
#                          y=ResidsSummStation_m1['WTD_MEAN'][select_iter].values.astype(float),
#                          name= m1,
#                          mode='markers',
#                          marker=dict(color =col1,
#                          size=mark_size,),
#                          ),
#                         row=1, col=1,
#                          )
# fig.add_trace(go.Scatter(x=ResidsSummStation_m2['STATION'][select_iter],
#                          y=ResidsSummStation_m2['WTD_MEAN'][select_iter].values.astype(float),
#                          name= m2,
#                          mode='markers',
#                          marker=dict(color =col2,
#                          size=mark_size,),
#                          ),
#                         row=1, col=1,
#                          )
# fig.add_trace(go.Scatter(x=ResidsSummStation_m3['STATION'][select_iter],
#                          y=ResidsSummStation_m3['WTD_MEAN'][select_iter].values.astype(float),
#                          name= m3,
#                          mode='markers',
#                          marker=dict(color =col3,
#                          size=mark_size,),
#                          ),
#                         row=1, col=1,
#                          )


# fig.add_trace(go.Scatter(x=ResidsSummStation_m1['STATION'][select_iter],
#                              y=ResidsSummStation_m1['WTD_RMS'][select_iter].values.astype(float),
#                          name= m1,
#                          mode='markers',
#                          marker=dict(color=col1,
#                          size=mark_size,),
#                          showlegend=False,
#                          ),
#                         row=2, col=1,
#                          )
# fig.add_trace(go.Scatter(x=ResidsSummStation_m2['STATION'][select_iter],
#                          y=ResidsSummStation_m2['WTD_RMS'][select_iter].values.astype(float),
#                          name= m2,
#                          mode='markers',
#                          marker=dict(color=col2,
#                          size=mark_size,),
#                          showlegend=False,
#                          ),
#                         row=2, col=1,
#                          )
# fig.add_trace(go.Scatter(x=ResidsSummStation_m3['STATION'][select_iter],
#                          y=ResidsSummStation_m3['WTD_RMS'][select_iter].values.astype(float),
#                          name= m3,
#                          mode='markers',
#                          marker=dict(color=col3,
#                          size=mark_size,),
#                          showlegend=False,
#                          ),
#                         row=2, col=1,
#                          )


# fig.update_layout(
#     title="Weighted Stats -- Residual Summary by Station",
# #     yaxis_title="Root mean Square",
# #     xaxis_title="Stations",
#     )

# fig.update_yaxes(title_text="mean residuals", row=1, col=1)
# fig.update_yaxes(title_text="rms", row=2, col=1)

# fig.update_xaxes(title_text="Stations", row=1, col=1)
# fig.update_xaxes(title_text="Stations", row=2, col=1)

# fig.update_layout(legend= {'itemsizing': 'constant'})
# fig.update_xaxes(tickangle=45)

# fig.update_yaxes(type="log", exponentformat= 'power')
# fig.update_layout(
#     autosize=False,
#     width=900,
#     height=700,
# )
# fig.update_layout(
#     font=dict(
# #             family="Courier New, monospace",
#         size=17,
#              ),)
# iplot(fig)

[ ]:

[ ]:

[ ]:

[ ]:

[ ]:

[ ]:

[ ]:

[ ]: